function [out_out_reordered,DAT_Trace,FractionatingMode_Out,bulkD] = YKG_wTRACE(inputfrommodel,elementsin,PINC,trace,Dmatrix)%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% This mfile is used to in conjunction with Kinzler and Grove 1992.% it is an adaptation of Yang et al. 1996 from FORTRAN to MATLAB by PM% Gregg 2007.  It calculates fractional crystallization paths given% an initial mantle composition.  This version of the code iterates% to look at compositions from 4 kbar to 0.001 kbar%% INPUT:    input - mantle composition.  must be 12 component and%               in the following order:%               SIO2, TIO2, AL2O3, FE2O3, FEO, MGO, CAO, NA2O, K2O,P2O5,%               CR2O3, MNO%           outfilenm - this is a character string for the output file%               which contains the oxide wt of the residual melt for each%               fractional xlization step% Outputs can vary as needed.% OUTPUT:   In the current version of the code outputs are all written to%           output files.  These files can vary based on what information%           the user wants outputed.%%% T. Gregg 2007%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%input = [46.58 0.95 17.66 0 9.31 11.90 11.43 1.81 0.37 0 0 0];%input = [48.62 0.94 18.31 0 7.92 10.09  9.93 3.82 0.37 0 0 0];%input = [49.37 0.89 17.49 0 7.98 10.28 10.46 3.31 0.23 0 0 0];% inputfrommodel% [SiO2,TiO2,Al2O3,Cr2O3,FeO,MgO,CaO,K2O,Na2O] wt%;%% Convert input elements into the correct oxide order from original Yang codeYKGelements = {'SIO2' 'TIO2' 'AL2O3' 'FE2O3' 'FEO' 'MGO' 'CAO' 'NA2O' 'K2O' 'P2O5' 'CR2O3' 'MNO'};input = ReorderColumnsOnly(inputfrommodel,elementsin,YKGelements);input(isnan(input))=0; input = input./sum(input)*100;%input = [inputfrommodel(1:3) 0 inputfrommodel(5:7) inputfrommodel(9) inputfrommodel(8) 0 inputfrommodel(4) 0];%%out = [];% KDOL = 0.3;% KDCPX = 0.23;% KDOL = 0.3;% KDCPX =.94*KDOL;%0.2726;%KDOL = 0.23;KDOL = 0.3;KDCPX = 0.23;%PINC = [1];%Main output file%file = ['runs-05.11.09/R2Du5Nu41e23_refinedT1325WH_YangOut_04kbar'];%Find Fe8 and Na8file2 = ['runs-05.11.09/Na8Fe8_out'];bulkD = [];kk = 0;for j = 1:1    %fid6 = fopen([file2(j,:),'.m'],'wt');    for i = 1:length(PINC)        kk = kk+1;        DAT = zeros(5,12);                                DAT(1,1:12) = input(j,:);                DAT_Trace(1,:) = trace;        FractionatingMode_Out(1,:) =[0 0 0 0 0 0];                %fid8 = fopen([file(kk,:),'.m'],'wt');        P = PINC(i);        %fprintf(fid8,['Out= [']);        %fprintf(fid6,'P(');fprintf(fid6,'%2.0f',i);fprintf(fid6,') = ');        %fprintf(fid6,'%8.4f',P);fprintf(fid6,';');                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        %  CRYSTAL FRACTIONATION OR ADDITION        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                CHECK1 = 0;        %CHECK1 = input(['\nENTER 0: FOR CRYSTAL FRACTIONATION \n',...        %     '      1: FOR OL-PL-AUG ADDITION. \n',...        %     '\n NOTE: THIS PROGRAM DOES NOT PERFORM OL OR OL-PL ADDITION. \n',...        %     'HUAI-JEN YANG HAS ANOTHER PROGRAM FOR DOING THESE. JUST ASK.\n']);        % (CHECK1.EQ.1) GOTO 500                %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        %  EVALUATING OLIV OR PLAG SATURATED        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                [QTZ,PLAG,OLIV,CPX,QQTZ,QOLIV,QCPX,PPLAG,POLIV,PCPX] = Tormey(DAT(1,:));        [ROQTZ,ROOLIV,ROPLAG,ROCPX,PROPLAG,PROOLIV,PROCPX,QROQTZ,...            QROOLIV,QROCPX] = OPALM_TRACE(DAT,P);        [XOPALM,YOPALM] = PLOT_TERN_TRACE(PROOLIV,PROCPX);        [XMODE,YMODE] = PLOT_TERN_TRACE(POLIV,PCPX);        DIFF1=XOPALM-XMODE;                test = 1;        if DIFF1 > 0            %%% 100            % Olivine fractionation            while (DIFF1 > 0)                [DAT(2,1),DAT(2,5),DAT(2,6)] = OLCOM_TRACE(DAT,KDOL);                                for J=1:12                    DAT(1,J)=(DAT(1,J)-DAT(2,J)*0.01)/0.99;                end                                FractionatingMode = [0 0 1 0 0 0];                D = FractionatingMode*Dmatrix';                                bulkD(:,test+1) = [D];                %D = Dmatrix*FractionatingMode';                DAT_Trace(test+1,:) = (DAT_Trace(test,:) - 0.01.*(DAT_Trace(test,:).*D))./0.99;                FractionatingMode_Out(test+1,:) =FractionatingMode;                                out = [out;DAT(1,1:12)];                %fprintf(fid8,'%8.4f',DAT(1,1:12));                %fprintf(fid8,' \n');                [QTZ,PLAG,OLIV,CPX,QQTZ,QOLIV,QCPX,PPLAG,POLIV,PCPX] = Tormey(DAT(1,:));                [ROQTZ,ROOLIV,ROPLAG,ROCPX,PROPLAG,PROOLIV,PROCPX,QROQTZ,...                    QROOLIV,QROCPX] = OPALM_TRACE(DAT,P);                [XOPALM,YOPALM] = PLOT_TERN_TRACE(PROOLIV,PROCPX);                [XMODE,YMODE] = PLOT_TERN_TRACE(POLIV,PCPX);                DIFF1=XOPALM-XMODE;                test = test+1;                            end                        %SMB: is this because it does an extra step?  should the previous            %step just be deleated?            %         for J=1:12            %             DAT(1,J)=(DAT(1,J)+DAT(2,J)*0.01)/1.01;            %         end            %         out = [out;DAT(1,1:12)];            %test            %fprintf(fid8,'\n');            %fprintf(fid8,'%8.4f',DAT(1,1:12));                    else            %%% 21            % Plagioclase fractionation                                    while(DIFF1 < 0)                                [DAT(3,1),DAT(3,3),DAT(3,7),DAT(3,8)] = OPAMPL_TRACE(DAT,P);                                for J=1:12                    DAT(1,J)=(DAT(1,J)-DAT(3,J)*0.01)/0.99;                end                out = [out;DAT(1,1:12)];                                                FractionatingMode = [0 0 0 1 0 0];                D = FractionatingMode*Dmatrix';                bulkD(:,test+1) = [D];                %D = Dmatrix*FractionatingMode';                DAT_Trace(test+1,:) = (DAT_Trace(test,:) - 0.01.*(DAT_Trace(test,:).*D))./0.99;                FractionatingMode_Out(test+1,:) =FractionatingMode;                                                %fprintf(fid8,'\n');                %fprintf(fid8,'%8.4f',DAT(1,1:12));                [QTZ,PLAG,OLIV,CPX,QQTZ,QOLIV,QCPX,PPLAG,POLIV,PCPX] = Tormey(DAT(1,:));                [ROQTZ,ROOLIV,ROPLAG,ROCPX,PROPLAG,PROOLIV,PROCPX,QROQTZ,...                    QROOLIV,QROCPX] = OPALM_TRACE(DAT,P);                [XOPALM,YOPALM] = PLOT_TERN_TRACE(PROOLIV,PROCPX);                [XMODE,YMODE] = PLOT_TERN_TRACE(POLIV,PCPX);                DIFF1=XOPALM-XMODE;                                test = test+1;            end                        %            %         for J=1:12            %             DAT(1,J)=(DAT(1,J)+DAT(3,J)*0.01)/1.01;            %         end        end                %% olivine and plag fractionation        A=100;        %%% 200        B = 1;        %test = 0;        while B > 0;            [DAT(2,1),DAT(2,5),DAT(2,6)] = OLCOM_TRACE(DAT,KDOL);            [DAT(3,1),DAT(3,3),DAT(3,7),DAT(3,8)] = OPAMPL_TRACE(DAT,P);                                    [ROQTZ,ROOLIV,ROPLAG,ROCPX,PROPLAG,PROOLIV,PROCPX,QROQTZ,...                QROOLIV,QROCPX] = OPALM_TRACE(DAT,P);                        [XOPALM,YOPALM] = PLOT_TERN_TRACE(PROOLIV,PROCPX);                        PLPROP=XOPALM/2*1.732;            OLPROP=1-PLPROP;                        PLPROP=PLPROP*2.7;            OLPROP=OLPROP*3.2;                        PLPROP=PLPROP/(PLPROP+OLPROP);            OLPROP=1-PLPROP;                        for J=1:12;                DAT(1,J)=(DAT(1,J)-(DAT(2,J)*OLPROP+DAT(3,J)*PLPROP)...                    *0.01)/0.99;            end                        out = [out;DAT(1,1:12)];                                    FractionatingMode = [0 0 OLPROP PLPROP 0 0];            D = FractionatingMode*Dmatrix';            bulkD(:,test+1) = [D];            DAT_Trace(test+1,:) = (DAT_Trace(test,:) - 0.01.*(DAT_Trace(test,:).*D))./0.99;            FractionatingMode_Out(test+1,:) =FractionatingMode;                                                %fprintf(fid8,'\n');            %fprintf(fid8,'%8.4f',DAT(1,1:12));                                              [QTZ,PLAG,OLIV,CPX,QQTZ,QOLIV,QCPX,PPLAG,POLIV,PCPX] = Tormey(DAT(1,:));                        [HJOLIV,HJPLAG,HJCPX,HJQTZ,QHJQTZ,QHJOLIV,QHJCPX,PHJPLAG,...                PHJOLIV,PHJCPX,MAL2O3,MCAO,MMGO,HJAL2O3,HJCAO,HJMGO] ...                = OPAM (DAT,P);                        [XQOPAM,YQOPAM] = PLOT_TERN_TRACE(QHJOLIV,QHJCPX);            [XMODE,YMODE] = PLOT_TERN_TRACE(QOLIV,QCPX);            DIFF2=sqrt((XQOPAM-XMODE)^2+(YQOPAM-YMODE)^2);            B=A-DIFF2;            if (B > 0);                A=DIFF2;            end            test = test+1;        end        %test        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%        % OLIVINE, PLAGIOCLASE AND AUGITE FRACTIONATION        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                        %SMB: is this because it does an extra step?  should the previous        %step just be deleated?        %     for J=1:12;        %         DAT(1,J)=(DAT(1,J)+(DAT(2,J)*OLPROP+DAT(3,J)*PLPROP)...        %             *0.01)/1.01;        %     end        %     for J=1:12;        %         DAT(1,J)=(DAT(1,J)+(DAT(2,J)*OLPROP+DAT(3,J)*PLPROP)...        %             *0.01)/1.01;        %     end                %%% 500        [QTZ,PLAG,OLIV,CPX,QQTZ,QOLIV,QCPX,PPLAG,POLIV,PCPX] = Tormey(DAT(1,:));                [HJOLIV,HJPLAG,HJCPX,HJQTZ,QHJQTZ,QHJOLIV,QHJCPX,PHJPLAG,PHJOLIV,...            PHJCPX,MAL2O3,MCAO,MMGO,HJAL2O3,HJCAO,HJMGO] = OPAM_TRACE(DAT,P);                [ROQTZ,ROOLIV,ROPLAG,ROCPX,PROPLAG,PROOLIV,PROCPX,...            QROQTZ,QROOLIV,QROCPX] = OPALM_TRACE(DAT,P);                [QCPXQTZ,QCPXOLIV,QCPXCPX,PCPXPLAG,PCPXOLIV,PCPXCPX,CPXOLIV,CPXPLAG,...            CPXCPX,CPXQTZ] = Tormey(DAT(4,:));                %test = 0;        % FORB = input(['FRACTIONATION = 1, ADDITION = 2 \n']);        FORB = 1;        STEPS = 10;        %        %     PROCPX        %     PROOLIV        %     PROPLAG                                OLPROP1 = 0.07;        PLPROP1 = 0.5;        CPXPROP1 = 0.43;                        OLPROP2 = 0.07;        PLPROP2 = 0.5;        CPXPROP2 = 0.43;                KDCPX1=DAT(1,6)/DAT(1,5)*DAT(4,5)/DAT(4,6);                        %%% 300        %while CONTROL <= 90 || (FORB == 1 & DAT(1,6) > 4.0 & DAT(1,6) < 9.0)        CONTROL=1;        while CONTROL <= 90 || ( DAT(1,6) > 4.0 && DAT(1,6) < 9.0)            [DAT(2,1),DAT(2,5),DAT(2,6)] = OLCOM_TRACE(DAT,KDOL);            [DAT(3,1),DAT(3,3),DAT(3,7),DAT(3,8)] = OPAMPL_TRACE(DAT,P);            [DAT(4,1),DAT(4,2),DAT(4,3),DAT(4,5),DAT(4,6),...                DAT(4,7),DAT(4,8),DAT(4,13)] = OPAMCPX_TRACE(DAT,P,KDCPX);                        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%            %  OUTPUT FRACTIONATION RESULT            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                        if CONTROL < STEPS;                OLPROP=OLPROP1;                PLPROP=PLPROP1;                CPXPROP=CPXPROP1;            else                OLPROP=OLPROP2;                PLPROP=PLPROP2;                CPXPROP=CPXPROP2;            end            if FORB == 1;                for J=1:12                    DAT(1,J)=(DAT(1,J)-(DAT(2,J)*OLPROP+DAT(3,J)*PLPROP+...                        DAT(4,J)*CPXPROP)*0.01)/0.99;                end            else                %             for J=1:12                %                 DAT(1,J)=(DAT(1,J)+(DAT(2,J)*OLPROP+DAT(3,J)*PLPROP+...                %                     DAT(4,J)*CPXPROP)*0.01)/1.01;                %             end            end                        out = [out;DAT(1,1:12)];                        FractionatingMode = [CPXPROP 0 OLPROP PLPROP 0 0];            D = FractionatingMode*Dmatrix';            bulkD(:,test+1) = [D];            %D = Dmatrix*FractionatingMode';            DAT_Trace(test+1,:) = (DAT_Trace(test,:) - 0.01.*(DAT_Trace(test,:).*D))./0.99;            FractionatingMode_Out(test+1,:) =FractionatingMode;                                    test=test+1;            %fprintf(fid8,' \n');            %fprintf(fid8,'%8.4f',DAT(1,1:12));                                                            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%            % PLOT RESULT            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%                        [QTZ,PLAG,OLIV,CPX,QQTZ,QOLIV,QCPX,PPLAG,POLIV,PCPX] = Tormey(DAT(1,:));            [HJOLIV,HJPLAG,HJCPX,HJQTZ,QHJQTZ,QHJOLIV,QHJCPX,PHJPLAG,PHJOLIV,...                PHJCPX,MAL2O3,MCAO,MMGO,HJAL2O3,HJCAO,HJMGO] = OPAM_TRACE(DAT,P);            [ROQTZ,ROOLIV,ROPLAG,ROCPX,PROPLAG,PROOLIV,PROCPX,...                QROQTZ,QROOLIV,QROCPX] = OPALM_TRACE(DAT,P);            [QCPXQTZ,QCPXOLIV,QCPXCPX,PCPXPLAG,PCPXOLIV,PCPXCPX,CPXOLIV,...                CPXPLAG,CPXCPX,CPXQTZ] = Tormey(DAT(4,:));            CONTROL=CONTROL+1;        end                           %         %         if any(inputfrommodel) > 0% %             out_out = [out(:,1:3) out(:,11) out(:,5:7) out(:,9) out(:,8) ((out(:,6)/40.3044))./(((out(:,6)./40.3044)+((out(:,5)./71.8464))))  ((out(:,6)/40.3044))./(((out(:,6)./40.3044)+((out(:,5)./71.8464))))...%                 (out(:,8)+out(:,9))./(out(:,8)+out(:,9)+out(:,7)) out(:,7)./out(:,3)];%             %             inputfrommodel=[inputfrommodel       ((inputfrommodel(:,6)/40.3044))./(((inputfrommodel(:,6)./40.3044)+((inputfrommodel(:,5)./71.8464))))...%                 (inputfrommodel(:,8)+inputfrommodel(:,9))./(inputfrommodel(:,8)+inputfrommodel(:,9)+inputfrommodel(:,7)) inputfrommodel(:,7)./inputfrommodel(:,3)];%             %             %             if size(inputfrommodel,2)==12;%                 inputfrommodel(13) = inputfrommodel(:,7)./inputfrommodel(:,3);%             end%             out_out = [inputfrommodel; out_out];%             %            %             allindices = [0:size(out_out,1)-1];%             %             Fout =1-((1-.01).^(allindices));%             PINCall = PINC.*ones(size(Fout));%             FractionatingMode_Out = [FractionatingMode_Out PINCall' Fout'];%         else%             %             out_out= NaN*ones(1,13);%             FractionatingMode_Out = zeros(1,8);%             DAT_Trace=zeros(1,45);%             %         end                                     [out_out_reordered,yikes] =  ReorderColumnsOnly(out,YKGelements,elementsin);            out_out_reordered = [inputfrommodel; out_out_reordered];            out_out_reordered(out_out_reordered==0)=NaN;            out_out_reordered = sumMgNumNorm_PT(out_out_reordered,[1:size(out_out_reordered,2)],[98.5 101.5]);                                                          allindices = [0:size(out_out_reordered,1)-1];            Fout =1-((1-.01).^(allindices));                        PINCall = PINC.*ones(size(Fout));            FractionatingMode_Out = [FractionatingMode_Out PINCall' Fout'];                                                                 % %     fprintf(fid6,'\n');        % %     fprintf(fid6,['Mg8',file(kk,:),' = [']);        % %     fprintf(fid6,'%8.4f',Mg8_int);fprintf(fid6,'];');        % %     fprintf(fid6,'\n');        % %     clear Mg8        % %     fprintf(fid8,'];');        % %    fclose(fid8);        %%% 1000 END            end    %fclose(fid6);        endendfunction [X,Y] = PLOT_TERN_TRACE (OL,CPX)% Subroutine to be used in conjunction with Yang et al. 1996.% This code has been adapted from FORTRAN by PM Gregg, 2007.%% THIS SUBROUTINE RECALCULATES DATA POINTS ON TERNARIES INTO X-Y COORDINATE% SUBROUTINE PLOT (OL,CPX,X,Y)% REAL OL,CPX,X,YX=(2-CPX/100-2*OL/100)/1.732;Y=CPX/100;endfunction [QTZ,PLAG,OLIV,CPX,QQTZ,QOLIV,QCPX,PPLAG,POLIV,PCPX] = Tormey(DAT)% Subroutine to be used in conjunction with Yang et al. 1996.% This code has been adapted from FORTRAN by PM Gregg, 2007.% It calculates tormey componentsMSIO2=DAT(1)/60.09;MTIO2=DAT(2)/79.9;MAL2O3=DAT(3)/102*2;MFE2O3=DAT(4)/159.7*2;MFEO=DAT(5)/71.85;MMGO=DAT(6)/40.3;MCAO=DAT(7)/56.08;MNA2O=DAT(8)/62*2;MK2O=DAT(9)/94.2*2;MP2O5=DAT(10)/141.9*2;MCR2O3=DAT(11)/152*2;SUM=MSIO2-MCAO-(2*(MNA2O+MK2O))+MTIO2+(2*MP2O5)+MFE2O3;QTZ=(MSIO2-(MMGO+MFEO)/2-(1.5*MCAO)-(MAL2O3/4));QTZ=(QTZ-(2.75*(MNA2O+MK2O)));QTZ=(QTZ+(MCR2O3/2)+(MTIO2/2)+(2.5*MP2O5))/SUM;PLAG=(MAL2O3+MNA2O-MK2O)/(2*SUM);OLIV=((MFEO+MMGO-MTIO2-MCAO)/2+(MAL2O3-MNA2O-MK2O-MCR2O3)...    /4+(0.83333*MP2O5))/SUM;CPX=(MCAO-(MAL2O3/2)+((MK2O+MNA2O)/2)-(1.66667*MP2O5))/SUM;OX=((MCR2O3/2)+MTIO2+MFE2O3)/SUM;AP=MP2O5/(3*SUM);OR=MK2O/SUM;QTZ=QTZ;PLAG=4*PLAG;OLIV=2*OLIV;CPX=3*CPX;QQTZ=100*QTZ/(QTZ+OLIV+CPX);QOLIV=100*OLIV/(QTZ+OLIV+CPX);QCPX=100*CPX/(QTZ+OLIV+CPX);PPLAG=100*PLAG/(PLAG+OLIV+CPX);POLIV=100*OLIV/(PLAG+OLIV+CPX);PCPX=100*CPX/(PLAG+OLIV+CPX);endfunction [ROQTZ,ROOLIV,ROPLAG,ROCPX,PROPLAG,PROOLIV,PROCPX,QROQTZ,...    QROOLIV,QROCPX] = OPALM_TRACE(DAT,P)% Subroutine to be used in conjunction with Yang et al. 1996.% This code has been adapted from FORTRAN by PM Gregg, 2007.%% SUBROUTINE OPALM (DAT,P,ROQTZ,ROOLIV,ROPLAG,ROCPX,PROPLAG,%                   PROOLIV,PROCPX,QROQTZ,QROOLIV,QROCPX)% DIMENSION DAT(51,12)% REAL DAT(51,12),MMGO,MFEO,MG,TIO2,NAK,P% REAL ROQTZ,ROOLIV,ROPLAG,ROCPX% REAL QROQTZ,QROOLIV,QROCPX,PROPLAG,PROOLIV,PROCPXMMGO=DAT(1,6)/40.3;MFEO=DAT(1,5)/71.85;MG=MMGO/(MMGO+MFEO);TIO2=DAT(1,2);NAK=(DAT(1,8)+DAT(1,9))/(DAT(1,8)+DAT(1,9)+DAT(1,7));ROQTZ=0.2-0.012*(P-0.001)+0.00*(1-MG)-0.252*NAK-0.00*TIO2;ROOLIV=0.12+0.008*(P-0.001)+0.229*(1-MG)-0.335*NAK-0.009*TIO2;ROPLAG=0.419+0.012*(P-0.001)-0.135*(1-MG)+0.691*NAK-0.00*TIO2;ROCPX=0.261-0.008*(P-0.001)-0.069*(1-MG)-0.104*NAK+0.012*TIO2;QROQTZ=100*ROQTZ/(ROQTZ+ROOLIV+ROCPX);QROOLIV=100*ROOLIV/(ROQTZ+ROOLIV+ROCPX);QROCPX=100*ROCPX/(ROQTZ+ROOLIV+ROCPX);PROPLAG=100*ROPLAG/(ROPLAG+ROOLIV+ROCPX);PROOLIV=100*ROOLIV/(ROPLAG+ROOLIV+ROCPX);PROCPX=100*ROCPX/(ROPLAG+ROOLIV+ROCPX);endfunction [PLSIO2,PLAL2O3,PLCAO,PLNA2O] = OPAMPL_TRACE(DAT,P)% Subroutine to be used in conjunction with Yang et al. 1996.% This code has been adapted from FORTRAN by PM Gregg, 2007.%% SUBROUTINE OPAMPL_TRACE(DAT,P,PLSIO2,PLAL2O3,PLCAO,PLNA2O)% DIMENSION DAT(51,12)% REAL DAT(51,12),MSIO2,MKALO2,MNAALO2,MCAAL2O4,MCAO,MMGO,MFEO% REAL MTIO2,MPO25,SUM,SUM1,A,MPLAN,MPLAB% REAL MPLCAO,MPLNA2O,MPLAL2O3,MPLSIO2,PLCAO,PLNA2O,PLSIO2,PLAL2O3MSIO2   = DAT(1,1)/60.09;MKALO2  = DAT(1,9)/94.2*2;MNAALO2 = DAT(1,8)/62*2;MCAAL2O4= (DAT(1,3)/102*2-MKALO2-MNAALO2)/2;MCAO = DAT(1,7)/56-MCAAL2O4;MMGO=DAT(1,6)/40.3;MFEO=DAT(1,5)/71.85;MTIO2=DAT(1,2)/79.9;MPO25=DAT(1,10)/141.9*2;SUM=MSIO2+MKALO2+MNAALO2+MCAAL2O4+MCAO+MMGO+MFEO+MTIO2+MPO25;MSIO2=MSIO2/SUM;MKALO2=MKALO2/SUM;MNAALO2=MNAALO2/SUM;MCAAL2O4=MCAAL2O4/SUM;A=MCAAL2O4/(MNAALO2*MSIO2)*exp(11.1068-0.0338*P-4.4719* ...    (1-MNAALO2)*(1-MNAALO2)-6.9707*(1-MKALO2)*(1-MKALO2));% MPLCAO=(KD*MLIQCAO/MLIQNA2O/2)/((KD*MLIQCAO/MLIQNA2O/2)+1)MPLAN=A/(1+A);MPLAB=1-MPLAN;MPLCAO=MPLAN;MPLNA2O=MPLAB/2;MPLAL2O3=(2*MPLAN+MPLAB)/2;MPLSIO2=2*MPLAN+3*MPLAB;PLCAO=MPLCAO*56.08;PLNA2O=MPLNA2O*62;PLAL2O3=MPLAL2O3*102;PLSIO2=MPLSIO2*60.09;SUM1=PLCAO+PLNA2O+PLSIO2+PLAL2O3;PLCAO=PLCAO/SUM1*100;PLNA2O=PLNA2O/SUM1*100;PLAL2O3=PLAL2O3/SUM1*100;PLSIO2=PLSIO2/SUM1*100;endfunction [HJOLIV,HJPLAG,HJCPX,HJQTZ,QHJQTZ,QHJOLIV,QHJCPX,PHJPLAG,...    PHJOLIV,PHJCPX,MAL2O3,MCAO,MMGO,HJAL2O3,HJCAO,HJMGO] = OPAM_TRACE(DAT,P)% Subroutine to be used in conjunction with Yang et al. 1996.% This code has been adapted from FORTRAN by PM Gregg, 2007.%% SUBROUTINE OPAM (DAT,P,HJOLIV,HJPLAG,HJCPX,HJQTZ,QHJQTZ,QHJOLIV,%                  QHJCPX,PHJPLAG,PHJOLIV,PHJCPX,MAL2O3,MCAO,MMGO,%                  HJAL2O3,HJCAO,HJMGO)% DIMENSION DAT(51,12)% REAL DAT(51,12),P% REAL HJQTZ,HJOLIV,HJPLAG,HJCPX,HJAL2O3,HJCAO,HJMGO% REAL QHJQTZ,QHJOLIV,QHJCPX,PHJPLAG,PHJOLIV,PHJCPX% REAL MSIO2,MTIO2,MAL2O3,MFE2O3,MFEO,MMGO,MCAO,MNA2O% REAL MK2O,MP2O5,MCR2O3,MTOTAL,OX,AP,OR% REAL SUM,QTZ,OLIV,CPX,PLAG,QQTZ,QOLIV,QCPX% REAL PPLAG,POLIV,PCPX% clear P% P = 0.001;MSIO2=DAT(1,1)/60.09;MTIO2=DAT(1,2)/79.9;MAL2O3=DAT(1,3)/102*2;MFE2O3=DAT(1,4)/159.7*2;MFEO=DAT(1,5)/71.85;MMGO=DAT(1,6)/40.3;MCAO=DAT(1,7)/56.08;MNA2O=DAT(1,8)/62*2;MK2O=DAT(1,9)/94.2*2;MP2O5=DAT(1,10)/141.9*2;MCR2O3=DAT(1,11)/152*2;MTOTAL=MSIO2+MTIO2+MAL2O3+MFE2O3+MFEO+MMGO+MCAO+MNA2O+MK2O+MP2O5+MCR2O3;MSIO2=MSIO2/MTOTAL;MTIO2=MTIO2/MTOTAL;MAL2O3=MAL2O3/MTOTAL;MFE2O3=MFE2O3/MTOTAL;MFEO=MFEO/MTOTAL;MMGO=MMGO/MTOTAL;MCAO=MCAO/MTOTAL;MNA2O=MNA2O/MTOTAL;MK2O=MK2O/MTOTAL;MP2O5=MP2O5/MTOTAL;%MCR2O3=CR2O3/MTOTAL;  < -- this line didn't work.  :P changed to next lineMCR2O3=MCR2O3/MTOTAL;HJAL2O3=0.236435+0.002182*P+0.109199*MNA2O+0.593295*MK2O-...    0.349718*MTIO2-0.299419*MFEO-0.130353*MSIO2;HJMGO=-0.276789+0.001139*P-0.542727*MNA2O-0.947294*MK2O-0.117029*MTIO2...    -0.489758*MFEO+2.085584*MSIO2-2.400115*MSIO2^2;HJCAO=1.132807-0.003388*P-0.569004*MNA2O-0.776108*MK2O-0.672285*MTIO2...    -0.213801*MFEO-3.355327*MSIO2+2.830219*MSIO2^2;SUM=MSIO2-HJCAO-(2*(MNA2O+MK2O))+MTIO2+(2*MP2O5)+MFE2O3;QTZ=(MSIO2-(HJMGO+MFEO)/2-(1.5*HJCAO)-(HJAL2O3/4));QTZ=(QTZ-(2.75*(MNA2O+MK2O)));QTZ=(QTZ+(MCR2O3/2)+(MTIO2/2)+(2.5*MP2O5))/SUM;PLAG=(HJAL2O3+MNA2O-MK2O)/(2*SUM);OLIV=((MFEO+HJMGO-MTIO2-HJCAO)/2+(HJAL2O3-MNA2O-MK2O-MCR2O3)/...    4+(0.83333*MP2O5))/SUM;CPX=(HJCAO-(HJAL2O3/2)+((MK2O+MNA2O)/2)-(1.66667*MP2O5))/SUM;OX=((MCR2O3/2)+MTIO2+MFE2O3)/SUM;AP=MP2O5/(3*SUM);OR=MK2O/SUM;HJQTZ=QTZ;HJPLAG=4*PLAG;HJOLIV=2*OLIV;HJCPX=3*CPX;QHJQTZ=100*HJQTZ/(HJQTZ+HJOLIV+HJCPX);QHJOLIV=100*HJOLIV/(HJQTZ+HJOLIV+HJCPX);QHJCPX=100*HJCPX/(HJQTZ+HJOLIV+HJCPX);PHJPLAG=100*HJPLAG/(HJPLAG+HJOLIV+HJCPX);PHJOLIV=100*HJOLIV/(HJPLAG+HJOLIV+HJCPX);PHJCPX=100*HJCPX/(HJPLAG+HJOLIV+HJCPX);endfunction [OLSIO2,OLFEO,OLMGO] = OLCOM_TRACE(DAT,KD)% Subroutine to be used in conjunction with Yang et al. 1996.% This code has been adapted from FORTRAN by PM Gregg, 2007.%% SUBROUTINE OLCOM_TRACE(DAT,KD,OLSIO2,OLFEO,OLMGO)% DIMENSION DAT(51,12)% 	REAL DAT(51,12),KD,MOLMGO,MOLFEO,OLMGO,OLFEO,OLSIO2,SUM% C	WRITE (9,50) (DAT(1,J),J=1,12)% 50      FORMAT (12F7.2)% moles MgO in OlMOLMGO=1/((KD*((DAT(1,5)/71.85)/(DAT(1,6)/40.3)))+1);%MOLMGO=1/((KD*DAT(1,5)/71.85/(DAT(1,6)/40.3))+1)% moles FeO in OlMOLFEO=1-MOLMGO;OLMGO=MOLMGO*40.3;OLFEO=MOLFEO*71.85;OLSIO2=0.5*60.09;SUM=OLMGO+OLFEO+OLSIO2;OLMGO=OLMGO/SUM*100;OLFEO=OLFEO/SUM*100;OLSIO2=OLSIO2/SUM*100;endfunction [CPXSIO2,CPXTIO2,CPXAL2O3,CPXFEO,CPXMGO,CPXCAO,CPXNA2O,TEMP] = ...    OPAMCPX_TRACE(DAT,P,KDCPX)% Subroutine to be used in conjunction with Yang et al. 1996.% This code has been adapted from FORTRAN by PM Gregg, 2007.%% SUBROUTINE OPAMCPX (DAT,P,KDCPX,CPXSIO2,CPXTIO2,CPXAL2O3,%                       CPXFEO,CPXMGO,CPXCAO,CPXNA2O,TEMP)% DIMENSION DAT(51,12)% REAL DAT(51,12),P,KDCPX,TEMP% REAL MSI,MTI,MAL,MFE2O3,MFE,MMG,MCA,MNA,MK,MP,MCR,MMN% REAL MCPXSI,MCPXTI,MCPXAL,MCPXALT,MCPXALO,MCPXNA% REAL MCPXMGFE,MCPXMG,MCPXFE,MCPXCA% REAL CPXSIO2,CPXTIO2,CPXAL2O3,CPXFEO,CPXMGO,CPXCAO,CPXNA2O% REAL MTOTAL,TOTALTEMP=0.4389*DAT(1,7)+4.14*DAT(1,8)-3.103*DAT(1,2);TEMP=TEMP+976;TEMP=TEMP+4.756*P;TEMP=TEMP+21.715*DAT(1,6);TEMP=TEMP+2.509*DAT(1,5);MSI=DAT(1,1)/60.09;MTI=DAT(1,2)/79.9;MAL=DAT(1,3)/102*2;MFE2O3=DAT(1,4)/159.7*2;MFE=DAT(1,5)/71.85;MMG=DAT(1,6)/40.3;MCA=DAT(1,7)/56.08;MNA=DAT(1,8)/62*2;MK=DAT(1,9)/94.2*2;MP=DAT(1,10)/141.9*2;MCR=DAT(1,11)/152*2;MMN=DAT(1,12)/70.94;MTOTAL=MSI+MTI+MAL+MFE2O3+MFE+MMG+MCA+MNA+MK+MP+MCR+MMN;MSI=MSI/MTOTAL*100;MTI=MTI/MTOTAL*100;MAL=MAL/MTOTAL*100;MFE2O3=MFE2O3/MTOTAL*100;MFE=MFE/MTOTAL*100;MMG=MMG/MTOTAL*100;MCA=MCA/MTOTAL*100;MNA=MNA/MTOTAL*100;MK=MK/MTOTAL*100;MP=MP/MTOTAL*100;MCR=MCR/MTOTAL*100;MMN=MMN/MTOTAL*100;MCPXCA=25.044+1.6237*MCA-0.83215*MFE-0.9531*MMG-0.6167*MAL+0.4457*(MNA+MK);% C	MCPXTI=-36.704+0.4661*MCA+0.3654*MFE+0.2590*MMG% C     +         +0.566*MAL+0.3825*(MNA+MK)+0.7905*MTI% C     +         +0.3020*MSI-0.034364*PMCPXTI=-3.6124+0.022523*MCA+0.0339*MFE+0.10776*MAL+0.0426*(MNA+MK)...    +0.3203*MTI-0.02169*P+0.001046*TEMP;% C        MCPXAL=-31.4919+0.5045*MCA+0.3309*MFE+1.6012*MAL% C     +         +1.2832*MTIMCPXNA=-2.2955+0.082*MNA+0.002214*TEMP-0.02024*MMG;MCPXMGFE=-11.051-1.0582*MCA+1.0893*MFE+1.1713*MMG+0.3848*MAL+0.5044*MSI;MCPXSI=116.082-0.31039*MSI-1.374*MTI-1.5822*MAL-0.5334*MFE-0.706*...    MCA-0.4177*MNA+0.2044*P-0.011125*TEMP;% C	MCPXALT=2*MCPXTI% C	MCPXALO=50-MCPXMGFE-MCPXCA-MCPXTI-MCPXNA% C        MCPXALT=50-MCPXSI% C	MCPXAL=MCPXALO+MCPXALTMCPXMG=((KDCPX*MFE/MMG)+1)/MCPXMGFE;MCPXMG=1/MCPXMG;MCPXFE=MCPXMGFE-MCPXMG;MCPXMG=MCPXMGFE*MCPXMG/(MCPXMG+MCPXFE);MCPXFE=MCPXMGFE*MCPXFE/(MCPXMG+MCPXFE);MCPXAL=98.2359-0.96422*MCPXSI-0.7958*MCPXTI-0.997656*MCPXMG-...    1.02176*MCPXFE-1.01482*MCPXCA-0.931973*MCPXNA;% C	MTOTAL=MCPXCA+MCPXTI+MCPXNA+MCPXMGFE+MCPXALO% C	WRITE (9,1001) MTOTAL,MCPXTI,MCPXAL% C1001    FORMAT (3F7.2)% C	MCPXSI=50-MCPXALT-MCPXTICPXSIO2=MCPXSI*60.09;CPXTIO2=MCPXTI*79.9;CPXAL2O3=MCPXAL*102/2;CPXFEO=MCPXFE*71.85;CPXMGO=MCPXMG*40.3;CPXCAO=MCPXCA*56.08;CPXNA2O=MCPXNA*62/2;TOTAL=CPXSIO2+CPXTIO2+CPXAL2O3+CPXFEO+CPXMGO+CPXCAO+CPXNA2O;CPXSIO2=CPXSIO2/TOTAL*100;CPXTIO2=CPXTIO2/TOTAL*100;CPXAL2O3=CPXAL2O3/TOTAL*100;CPXFEO=CPXFEO/TOTAL*100;CPXMGO=CPXMGO/TOTAL*100;CPXCAO=CPXCAO/TOTAL*100;CPXNA2O=CPXNA2O/TOTAL*100;end